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Abstract. 

This article discusses the existence of solitary electromagnetic waves trapped 
in a self-generated Langmuir wave and embedded in an infinitely long circularly 
polarized electromagnetic wave propagating through a plasma. From the 
mathematical point of view they are exact solutions of the 1-dimensional 
relativistic cold fluid plasma model with nonvanishing boundary conditions. 
Under the assumption of traveling wave solutions with velocity V and vector 
potential frequency uj, the fluid model is reduced to a Hamiltonian system. The 
solitary waves are homoclinic (grey solitons) or heteroclinic (dark solitons) orbits 
to fixed points. By using a dynamical systems description of the Hamiltonian 
system and a spectral method, we identify a great variety of solitary waves, 
including asymmetric ones, discuss their disappearance for certain parameter 
values, and classify them according to: (i) grey or dark character, (ii) the number 
of humps of the vector potential envelope and (iii) their symmetries. The solutions 
come in continuous families in the parametric V — uj plane and extend up to 
velocities that approach the speed of light. The stability of certain types of 
grey solitary waves is investigated with the aid of particle-in-cell simulations 
that demonstrate their propagation for a few tens of the inverse of the plasma 
frequency. 



PACS numbers: 52.35.Sb, 52.38.Kd 



1. Introduction 

The excitation of long-lived solitary waves during the interaction of high- intensity laser 
pulses with plasmas is a topic with applications and of theoretical interest. As multi- 
dimensional particle in cell (PIC) simulations have shown, these waves form behind 
the laser pulse and they consist of electron density depressions with a trapped intense 
electromagnetic field oscillating at a frequency well below the laser frequency [l]-[5j. 
Solitary waves can propagate towards the plasma-vacuum interface where the stored 
electromagnetic energy is radiated away in the form of low-frequency electromagnetic 
bursts [2j, a process recently detected in the laboratory Bright spots in optical 
plasma images with the same polarization as the laser pulse have been attributed to 
the formation of such solitons [7j. These waves can evolve to a state named postsoliton, 
that has also been observed in the laboratory with proton imaging techniques [8 -11 



The propagation of electromagnetic solitary waves has been intensively 



investigated within the cold, relativistic, one-dimensional fluid approximation 12 -22 



Restricting attention to circularly polarized travelling wave solutions, with velocity V 
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and frequency uj, yields a pair of second order differential equations that governs the 
dynamics of the electrostatic potential (j) and the amplitude of the vector potential a. 
Hence, the system describes a Langmuir and an electromagnetic wave coupled by the 
nonlinear terms arising from the perturbation of the density and the relativistic mass. 
It admits solitary waves with vanishing (VBC) and nonvanishing (NVBC) boundary 
conditions and, even though the governing equations are not completely integrable, 
they are commonly referred to as solitons. 

For VBC, a — ?> and as a: ~> ±oo, a soliton is interpreted as a light 
wave which is trapped in a self-generated plasma wave 13 . They have been classified 



according to the number of zeros of the vector potential profile, p, see Ref. 17 . They 
are commonly referred to as bright solitons and their stability has been analyzed 
too [4{ [l8{[20|[23}|28] . In 17 , the existence of such solutions is discussed in the 
parametric V — ui plane and p = 0,1, 2... families of solitons are identified. Some of 
these families end at certain velocity values where the ion density profile shows a cusp 
at the center of the soliton. The soliton breaking has been proposed as a mechanism for 
particle acceleration in high-intensity laser plasma interaction l7 and it was observed 
in PIC simulations for solitons with V — and overcritical amplitude 15 . Solitons 
with VBC in warm 



20 and magnetized plasmas 16 22 have been studied too. 



On the other hand, three different solutions with NVBC are possible 19 21]: (i) 
grey solitons (a — > oq and — > as a; — > ±oo), (ii) dark solitons (a — >■ ±ao and </>—>■ 
as a; — >■ ±oo) and (iii) shock waves (a — >■ oq and ^ — >■ as a — — oo and a — >■ 
and — >■ (/)o as X — >■ -l-oo). The asymptotic values oq and 0o are related to the two 
parameters V and ut. The branch of shock waves splits the V — ui plane in two different 
regions where either dark or grey solitons exist. Branches of solutions are found and 



they break down at increasing oq due to divergence of the electron density 19 21 



As we will see, the grey and dark solitons with NVBC can be interpreted as 
a localized modulation in a long circularly polarized electromagnetic wave coupled 
with a plasma wave. Even though the circularly polarized wave is susceptible to 
the relativistic Raman and the modulational instabilities 29 , the analysis of these 



solitary structures is fully justified. First, the present work extends the discussion 
about solitary waves with NVBC, that until now was limited to the narrow velocity 
range, < V/c < 0.051 19 21 . Second, the solitary waves could play a role in 
processes that are faster than the inverse of the parametric instability growth rate. 
We also point out that any damping mechanism, not included in the present analysis 
for simplicity, could reduce or suppress these instabilities. 

From the mathematical point of view, the problem of finding solitons reduces, 
through the traveling wave ansatz, to that of finding orbits of an associated 
Hamiltonian system for a and 4> consistent with the boundary conditions of the fiuid 
model for each type of soliton. The Hamiltonian system is four dimensional, time- 
reversible and autonomous. It admits four fixed points that, for brevity, will be 
denoted by the letter Q and the value of </> and a inside brackets. Only three of them, 

= (0,±ao) and Qi — {(f)o,0), play a role when discussing the existence of solitary 
waves. Orbits that are asymptotic to an equilibrium as a: — >■ ±oo are called homoclinic 
connections, while orbits that connect two distinct equilibria in the limits x — >■ ±oo are 
termed heteroclinic connections (see 30 for a review of homoclinic orbits in reversible 
systems). Therefore, grey solitons, dark solitons and shock waves correspond to the 
connecting orbits Qq — )■ Qo^j Qq Qo and Qq Qi respectively. 

Existence and robustness under variations in V and w of such connecting orbits 
can be directly inferred by simple geometrical arguments from the theory of dynamical 
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systems, with a critical role played by the properties of the equilibria of the system 



and in particular their linear stability. Previous works 19 21 studied the parametric 
domain where is a saddle-center, with eigenvalues A1.2 = ia, A3. 4 = ±i/3, and 
they found solutions in the velocity range < V/c < 0.051. 

Here, we extend the analysis to velocities up to V/c = 1, and specifically in regions 
of the parametric domain where is a saddle-focus, i.e. its eigenvalues are of the 



form ±(a ± iP). It is well known 31 that homoclinic connections of a saddle- focus 
in Hamiltonian systems do generally exist and are robust against small perturbations 
of the Hamiltonian. Therefore one expects them to appear in continuous families in 
the V — Lu plane. Moreover, existence of such a homoclinic orbit, for instance one 
that corresponds to a one-hump soliton, implies existence of infinitely many of them, 
corresponding to solitons with a different number of humps 31 . Taking into account 
these insights from the theory of dynamical systems, we numerically identify new 
families of single and multi-hump solitons of the grey and dark varieties that can be 



interpreted as counterparts to the branches reported for VBC 17 . We also report, 
for the first time, asymmetric solitons in this system and explain the disappearance 
of certain families of solutions as parameters are varied. 

The organization of the paper is as follows. In section |2] the system of equations 
that governs the dynamics of the solitons is revisited, while its fixed points and their 
possible connecting orbits are discussed in section [3j Section |4] describes the spectral 
algorithm 32 that we have used to find homoclinic and heteroclinic orbits, as well as 
some optimizations carried out for our specific problem. The numerically computed 
families of dark and grey, symmetric and asymmetric solitons are presented in section 
|5] In section |6] the stability of certain types of grey solitons is investigated with the 
aid of particle-in-cell (PIC) simulations. The conclusions are summarize in section [t] 
where the similarities and differences with the VBC waves are stressed. 



2. Dynamical equations 

This section briefiy summarizes the theory of one-dimensional circularly polarized 



solitons (see 12 17 33 for a detailed discussion). The plasma is assumed to be cold 
and composed of electrons and ions, denoted by the subscript a — e,i respectively. 
It is convenient to use length, time, velocity, momentum, vector and scalar potential, 
and density normalized over c/ujpe, c, niaC, mgC^/e and ng respectively. Here 

no, rria and ujpe — (47rrtoe^/me)^/^ are the unperturbed density, the rest mass and the 
electron plasma frequency. Using this notation the Maxwell (in the Coulomb gauge) 
and plasma equations read 

d^A d 

AA - - — V0 = rigVe - n^v^ (la) 
A(j) = rie — Ui (16) 

^+V-(n„v„)=0 (Ic) 
dV 

^ - v„ X (V X P„) = -V(e„0-H7a) (Id) 

where A and (j) are the vector and scalar potential, Pq, = Pa+eaA, 7q = (1 + |PqP)^^^ 
and Pq, and = TPa/la a-re the kinetic momentum and the fluid velocity respectively. 
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For convenience the dimensionless parameter 
{qe — ~P' and qi — e are species charges). 



{qaTrie) / (enia) has been introduced 



Taking dy = dz = 0, the Coulomb Gauge immediately gives A — A± where _L 



denotes the direction perpendicular to x. The transverse component of (Id) yields 

P_La = 0. We also take all the variables to be functions of ^ = (x — Vt)/\/ 1 — V^ and 
assume a circularly polarized vector potential 

Ay + iAz^a{C)e'^^& (2) 

Here, V is the group velocity of the solitary wave, k is the wavevector and uj is the 
frequency in a frame moving with the solitary wave. With the above assumptions, 



equation (Ic) and the longitudinal component of (Id) can be integrated. Imposing the 



boundary conditions a = ±ao, (J) — 0, — 1 and p^a = as x — > — oo, the kinetic 
momentum, energy and the density of each species are functions of just the potentials 
(f) and a. For instance, the densities and the 7 factors are given by 

(1 - V^)r^ 



na{(t), a) 



(3) 



7a 



•0a - Vrg 

i-y2 



(4) 



with ipa = To 



r„ ^ [02 _ (1 „ v^){l + eW)] and F^ = (1 + eWoY/^ 
For brevity we write — )■ e and for the numerical calculation we set e = 1/1836. 



Substituting in (la) and (16) yields 19 21 



V 



1 



V 



(5 a) 



(56) 



and V — uj/k, where the prime denotes derivative with respect to ^. We remark 
that, for solutions described by equation ([2]), the group velocity of the solitary wave V 
should be equal to the phase velocity uj/k. Solutions with different velocities require 
the addition of a relative phase in equation ^ 21 



System (5a)-(56) describes the dynamics of localized electromagnetic modulations 
trapped by a self-generated Langmuir wave embedded in a infinitely long 
electromagnetic wave. Introducing the momenta Pa = {1 — V^)a' and — — 0', 
it can be written as a fourth order Hamiltonian system with Hamiltonian 



H{a,Pa,(l),P^) 



Pr, 



1 2 2 



1 9 



V 



(«,0) 



(6) 



34 



We also note that the conditions > imply the following restriction on 

v/(l - V^){1 + a2) - Te < <^ < J (r, - v/(l - V2)(l + 62a2)) (7) 
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3. Fixed points and connecting orbits 

3.1. Fixed points and their stability 



The fixed points are obtained by setting the righthand side of (5a)- (56) equal to zero 
and, for brevity, they will be denoted by Q{(t>, a) (assuming 0' = a' = as ^ — ?► ±oo). 
We note that only orbits connecting the fixed points = (0, ±ao) are consistent with 
the previously imposed boundary conditions. Hence, Qq and their stability properties 
play an important role when discussing the presence of solitary waves in the full fluid 
system. 

3.1.1. Fixed points Qq . One readily verifies that ~ (0,±ao) is a fixed point of 
system (5a)-(56) if the following dispersion relation is satisfied 

1 e 



(8) 



Its physical interpretation is evident when we express (p|) as a function of the 



wavevector k^p and the frequency in the laboratory frame ujlf 17 21 

1 e 



f^LF 



(9) 



that is the normalized dispersion relation of a pure transverse electromagnetic wave 
with relativistic effects and the ion motion. Note that equation ^ has a solution 
within the range < < 1 + e. 

The stability of is determined by the eigenvalues of the Jacobian of system 



(5a 



u} plane there exist 3 different regions of stability; for any given 
ui, Qf is a saddle-center if F < Vsc, a center if Vsc < V <Vsf and a _saddle-focus 



56|). In the 

_ ^" is a saddle-c( 

if V > VsF- The velocities Vsc and Vsf are (see Appendix [t] and figure [I]) 
Vsc = + 



(r? + eri)2 



[2ao(i-62)rer.]" 

-[r?(ri + a2) + eri(r? + e2ag)]^ 



Vsf 

Note that Vsf 



1 - 



2ao(i-62)rer, 



r?(rn-ag)-f er3(r? + e 



Vsc ^ e/(l 



Vsc 1 as ao 
e + e^) as flo (w^ ^ 1 + e) 



oo Iw 



0) and Vsf 



(10) 

(11) 
and 



3.1.2. Fixed point Qi. System (5a)-(56) also admits the fixed point Qi — ((Foi — 
roe)/(l -I- e),0). Since orbits must connect Qi to Qq , we need to enforce the same 
value of the invariant (equation (|6|) at both fixed points. This yields 

21-1/2 



Vs = 



1 - 



^0 



2r,r, 



l + e 



^0 



2rj,: 



(12) 



which gives an explicit relation between the velocity of the wave and the asymptotic 
value ao, plotted in figure [l] As ao H. 1 e) one has Vs ^ e/{l — e + e^). 
Appendix [t] shows that Qi is a saddle-center for velocities given by ( [121 ) ^^^^ ™ 
principle, shock waves are possible in the range e/{l — e + e^) < V < 1. Shock waves 
belong to the parametric domain where Qq is a saddle-center, see figure [ij and will 
not be studied in the present work (see 19 and [2l]). 
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3.1.3. Fixed point Q^. The system has an additional fixed point Q2 for velocities V 
greater than 

_ r, + 6r, ^(r, + eV^Y - (i + ey 

(see appendix [?]) . As Q J and Q2 must share the same value of H for a connection to 
exist, a condition similar to ( [l2| ) can be obtained, which cannot be fullfilled for any 
value of the parameters and therefore heteroclinic connections involving Q2 and Qq 
are not possible (see Appendix [t]) . 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



V 



Figure 1. Characteristic velocity curves and Qq stability regions in the V — u) 
plane. The velocities Vsci'^) a^nd Vsf{'^) (solid lines) split the plane in three 
regions where has different stability. Along the curve Vsi^jj) (dashed line) 
shock waves are possible. The inset shows a detail close to V ~ 1 where both 
VsF and Vs approach 1 as oj — > 0. 



3.2. Connecting orbits 

We now turn to the discussion of the conditions under which heteroclinic and 



homoclinic connections, and therefore solitons of the fluid system {la]-{ld] are 
expected to exist. We will need to recall well known facts from the theory of 
dynamical systems (see [30| , for example) and introduce the notion of the stable 
(unstable) manifold {W^) of a fixed point Qi, as the set of forward (backward) in 
^ trajectories that terminate at Q^. A heteroclinic (homoclinic) connection from 
to (to itself) lies on the intersection of the unstable manifold of Q' and the stable 
manifold of (Q*). 

Stable and unstable manifolds are complicated objects with intriguing structure 



and their visualization can become a formidable task (see 35 for a review of methods) 
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Here it will be enough to consider the dimensionality of stable and unstable manifolds, 
as this dictates whether or not they are in general expected to intersect. We note that 
dynamics of the 4-dimensional system (5a|-(56) are constrained to a 3-dimensional 



energy-manifold by conservation of the invariant H, equation ([6|. The intersection of 
stable and unstable manifolds has to take place within this energy-manifold. 

We first discuss homoclinic connections involving Q^. When Qq is a saddle- 
center, both its stable and unstable manifolds are 1-dimensional and are, in general, 
not expected to intersect in 3-dimensional space. If they do intersect, they necessarily 
have to be the same from uniqueness of solutions. A homoclinic connection is therefore 
not generic, in the sense that it requires the condition Wq — Wq which can only 
be fullfilled for specific V and uj values. Hence the homoclinic connections in the 
saddle-center case are expected to occur in branches in the V — uj plane. However, the 
existence of a continuous spectrum of single-hump solitons when is a saddle-center 
has been suggested for VBC 



18 



19 21 boundary conditions. Unless 



and NVBC 

further restrictions are imposed, such a continuous spectrum has to be interpreted 
as a numerical artifact, only valid within the accuracy of long-time integration of 
system (5a|-(56|. Nevertheless, singular cases could arise. For instance, at the V — 



case where the dynamics is 2-dimensional, it can be proved analytically that standing 
solitary wave solutions exist within a continuous uj range ^15^ . 

In the saddle-focus case on the other hand, Wq and Wq are 2-dimensional and 
are in general expected to intersect transversally along a 1-dimensional curve in the 
3-dimensional energy-manifold. Therefore we expect homoclinic solutions to exist 
generically in the V—uj plane. Furthermore it is well known that, for given parameters, 
if one such transverse intersection exists then there exist infinitely many, which have 
been shown to form a local, complete Horseshoe structure 31 . In our case this implies 
that for any one-hump soliton we can find an associated family of multi-hump solitons 
for any given V and w [36) . 

The symmetry properties of (5a|-(56| help us deduce more properties of the 
connecting orbits. Equations (5a)-(56| arc invariant under the reflection symmetry, 
a — ^ —a, a' — >■ —a' . As a result, for any homoclinic Qq — >■ Qq connection there exists 
an identically shaped Qq — > Qq one. Note that, since we impose inhomogeneous 
boundary conditions, we cannot have reflection invariant connections. Moreover, 
reflection symmetry allows us to carry many of the results available for homoclinic 
orbits over to heteroclinic orbits involving Qq and Qq (dark solitons); the two 
equilibria can be considered as a single equilibrium in a reduced system in which 
a is identified with —a and in which heteroclinic connections become homoclinic 36 



Therefore we can expect heteroclinic orbits in the saddle-focus parametric regime to 
occur generically. Furthermore, the Qq — > Qq connection has an identically shaped 
Qq — > Qq counterpart. 



System (5a)-(5&l is also time-reversible, that is invariant under simultaneous 
change of sign of ^ and the generalized momenta P^, Pa- As a result, homoclinic orbits 
can either be symmetric (self-dual, (/)(^) = (/>(— C), a(0 — oi' come in pairs of 

asymmetric orbits related by time-reversal. Similarly, heteroclinic connections Qq — >■ 
Qq have to be either antisymmetric functions of ^ (0(0 — 0(— C)) «(0 = 
or come in asymmetric pairs, related by ^-reversal. Note that the antisymmetric 
heteroclinic solitons are associated with the combined action of ^reversal and 
reflection. Asymmetric homoclinic and heteroclinic orbits in ^-reversible Hamiltonian 
systems are well studied ^30] , but to our knowledge they appear here for the first time 
in the context of relativistic solitons, see section [S] The addition of a perturbation 
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that breaks the conserved quantity H but preserves the reversibihty would destroy 
asymmetric solutions 30 . 

For later use, we introduce the so-called symmetric section S : a' = (j)' = 0, which 
is in a sense a symmetry hyperplane of the system, as any point in it is left invariant 
under ^-reversal. Symmetric homoclinic orbits have to intersect S 



30 . As we will see 



in section [5] symmetric homoclinic orbits can disappear through a mechanism referred 
to as coalescence 30 , as a result of which the stable and unstable manifolds fail to 



intersect on the symmetric section S. 



4. Numerical algorithm 



Due to time reversibility, the computation of symmetric homoclinic or antisymmetric 



heteroclinic orbits of (5al 



56| involves the solution of a boundary value problem on 
Within the regime where Qq is a saddle-center and 
previous works truncate this interval and 



the semi- infinite interval (—00, 0) 
the unstable manifold is one-dimensional 
consider the linearized dynamics close to the fixed points, solving for the parameter 
values in the V—u) plane in which the boundary value problem is satisfied 12|l7|19pr 



For the saddle-focus case, where the local unstable manifold is a 2-dimensional plane, 
one would also need to solve for the initial condition in this plane, for instance by 



parametrizing it by a polar angle 37 



Here we implement the rational spectral collocation algorithm of Ref. 32 that 
avoids both the truncation of the domain and the introduction of an additional 
parameter. Let us consider the system u" — f(u) with u = (a, 0)'^ and the 
components of f given by the right-hand sides of (5a)- (5&I. Assuming a fast enough 



±00 38 , the 



approach of the solutions to the asymptotic values u u± as 
variables can be expanded as a sum of orthogonal rational functions 

Af+l 

= ^^fc'^o^ [kcot'HO] , 1 = 1,2 (14) 

k=0 

The basis functions suggest the following choice for the M collocation points 



~ cot 



1 < j < Af 



(15) 



that are complemented by the two following collocation points ^0 = +00 and 
£,M+i = —00. The coefficients Cik are given by 

n M+i 



cos I I Q<k< M +1 (16) 



with 



2 if m = or TO = A/ -I- 1 and = 1 if 1 < w < M . 



u" = /(u) yield 2M nonlinear algebraic equations 



Computing the second derivative from ( 14 ) and substituting the results in 

i = l,2 (17) 



M+l 



DjkC^k + h [u(e,)] =0 1 < J < M, 



fc=0 



where 
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We also have 4 boundary conditions 

M+l 

C,k = u+ (19a) 



fc=0 



J2 = u- (196) 

fc=0 

with 1 < i < 2. 

Since the system is autonomous, every translation on ^ of its solution is also a 
solution. This indeterminacy can be removed by adding the phase condition 

<u'(e)-ii'(o,u"(e)>=o (20) 

where < •,• > denotes the P inner product and u the previous orbit on a branch. 
Solutions occuring at branches tu = ijj{V) need this extra condition to adequately vary 



the parameter. This would be the case of the branches found for V < Vsc 19 21 



However, as we will see, for V > Vsf the solutions exist continuously in the V 



plane. In this case equations (17l and (19a)-(19&) constitute a set of 2(M+2) equations 



for the 2{AI + 2) coefficients Cik that can be solved for fixed values of w and V. We 
used a Newton-Raphson algorithm with its Jacobian calculated analytically to speed 
up the convergence. 



Since system (5a)-(5&) does not include the first derivatives of the variables {(f), a) 



on the right-hand side, we have directly computed the second derivative. This is a 



difference from the general scheme given in reference 32 that allows us to divide 
by two the number of coefhcients. Note also that for asymmetric solutions all the 
coefhcients must be calculated whereas in the case of symmetric or antisymmetric 
solutions just one half of the coefficients are needed (for instance symmetric grey 
solitons are even solutions and Ci^2k+i — 0). 

The variable ^ can be stretched according to ^ — ?> L^. Here L is a scaling factor 
that can be used to optimize accuracy of the solution. Although some strategies 
have been discussed to find a proper value of L [39| , in practice it is chosen by 
experimentation with different L for a given value of M (as suggested in reference [40]). 
We found L = 6 to be an adequate value for our numerical computations and a 
number of collocation points M ^ 300. The initial guess for the Newton-Raphson 
algorithm can be obtained from the intersection of stable and unstable manifolds, 
as discussed in section [Sj For asymmetric solutions it can also be built up from 
pieces of symmetric solutions of different parity, smoothed by dropping the high 
order coefficients. Similarly, guesses for multi-hump solutions can be built up from 
single hump solutions of appropriate shape. Once we have a good initial guess, the 
spectral algorithm gives an accurate solution that can be used as initial guess for other 
parameter values. 

We finally remark that the accuracy of the solutions can be checked by testing the 
spectral convergence: a plot of the logarithm of the coeficients Cik versus the order 
k should decay linearly. This is an advantage as compared to the truncation of the 
interval methods where the error is controlled by the initial distance to the fixed point 
which is fixed arbitrarily. Such a criterion allows to separate true from false solutions, 
specially when the fixed point has eigenvalues with small real part (in absolute value) 
and the solutions approach slowly to their asymptotic values. 
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5. Numerical results 

5.1. Families of solitary waves 

Figure [2] depicts our numerically computed grey soliton solutions, belonging to a 
certain family. In panel [2ja) we also plot the curves Vsc and Vsf- The panels (d)- 
(f) display a representative solution at the parameter values V = 0.8 and uj — 0.11, 
just above the frequency uj ^ 0.10545 where these solutions disappear for V = 0.8. 
The densities and the 7 factors are computed from equations ([3| and Q respectively 
(note the different scales for 7e and 7^). We point out that for this particular family 
of solutions the ion density and the potential exhibit one hump at the center 
of the solution. The physical meaning of the solitary waves with NVBC can be 
undestood with the aid of figure [3] where we plotted electric field and electron momenta 
components. The longitudinal electric field and momentum p^e vanish outside the 
solitary wave that, as the plots of Ey and Pyf, show, can be intepreted as a localized 
modulation of an infinitely long circularly polarized electromagnetic wave {E^ and pze 
are not presented). 

Even though our fluid description of the plasma cannot take into account discrete 
particle effects like acceleration or heating, it is interesting to compute the maximum 
value of the electron and ion kinetic energy within the solitary wave 

£;™'^- = max[mc,c\j^ - 1)] (21) 

For the previously discussed family, these quantities are plotted in figure [2jb) and[2](c) 
(z-axis is in a logarithmic scale). The electron fiuid can reach energies of the order of 
hundreds of MeV and ions several tens of MeV. We also note that, for a fixed value 
of UJ, the largest energy is reached at the lowest admissible V value. 

Figure |4] corresponds to a different family of grey solitons. In this case both 
potentials present a minimum at the center of the soliton and the electron density 
has one central peak. The panel |4]ja) shows the existence domain found by randomly 
varying the parameters V and uj (as we also did for figure [2] ) . However, by fixing 
one of the parameters, using a small step for the other and initializing the algorithm 
with the solution obtained at the previous iteration, we were able to find solitary 
waves of the same type for parameter values outside the region indicated in |4|^a) (see 
figure |4]jd)-|4]^f)). As the frequency decreases, the potential develops a cusp shape at 
the center of the wave and the algorithm requires an initial condition closer to the 
real solution to ensure convergence. Hence, the boundary exhibited in figure Qa) is 
a numerical artifact as we will also confirm in section [5?2| We remark that potential 
profiles with a cusp shape have also been reported in solitary waves with VBC for 



parameter values close to the wavebreaking 17 . 

A representative family of dark solitons, or heteroclinic connections — Qq , is 
shown in figure [5j Note that its existence domain encompasses that of the grey solitons 
of figure[2] Panels [5]Jd)-[5]Jf) display a particular solution of this family of dark solitons 
at the parameter values ut — 0.08 and V = 0.8. The potential 4> has one central hump 
and, as opposed to the grey solitons, the vector potential is an antisymmetric function. 
The peak electron and ion kinetic energies are of the order of hundreds and tens of 
MeV, respectively. 

These three families of solitary waves provide some of the simplest examples of 



solitary solutions admitted by the system (5a)-(56). However, due to the fact that the 



fixed point is a saddle-focus and the system is Hamiltonian and reversible, multi-hump 
and asymmetric solutions can exist too. A few examples of grey multi-hump solutions 



Relativistic solitary waves in plasmas 



11 




Figure 2. (Color online) A family of one-hump symmetric grey solitons. 



with Lu = V = 0.8 are displayed in figure [6]J a)- [6]Jd). Panels |6je) and[6]^f) correspond 
to the density and 7 factor of the solution with 5 humps. The potential (f> exhibits 
one central hump, while the vector potential has multiple humps. These solutions 
are characterized by a cavity with a depression of the electron density in which an 
electromagnetic wave can be trapped. On the other hand, grey and dark asymmetric 
solitons are shown in figure [T^a) and[7jb), respectively For completeness. Figure [Tj^c) 
also shows a multi-hump dark soliton (heteroclinic connection). Solutions in panels 
[rja) andjrjc) are plotted in panel [7]jd), where a projection on the phase space (j) — a 
together with the fixed points Qq is shown. Note that the asymmetric grey soliton 
connects Qq — Qq, although it passes close to Qq (see the inset). 



5.2. Coalescence and disappearance of solitary waves in parameter space 

A striking feature of the families of solitary waves we identified is that they fill in 
regions with a well defined, within our numerical precision, boundary in the V — uj 
plane. As we have seen in section [3] homoclinic and heteroclinic orbits lie on the 
interesection of stable and unstable manifolds of fixed points. This suggests that 
certain families of orbits can cease to exist when the corresponding manifolds fail to 
intersect in the neighborhood of the solitary solution. The mechanism responsible for 
loss of intersection in our case is known as coalescence [30[|41[[42] and we now turn to 
its detailed description, through a numerical experiment involving the disapperance 
of a member of the family of grey solitons shown in figure [2] 

An approximation to the unstable manifold of can be computed by integrating 
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equations (5a|-(56| with initial condition 

[(j) (p a a] ^ [0 ao 0] + ctvi 



(22) 
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Figure 5. (Color online) A family of one-hump antisymmetric dark solitons 



where vi is the real part of one of the unstable eigenvectors of Qq , associated with an 
eigenvalue A = Xr + iXi, Xr > 0. Here cr is a small parameter that controls the position 
of the initial condition on the plane spanned by the unstable eigenvectors. Computing 
one thousand outward spiraling orbits with a — 0.001 xexp[(27rAr)/(jAi)], j ~ 1...1000 
then results in an approximation of the unstable manifold. The stable manifold was 
computed similarly, sprinkling initial conditions along the stable eigendirection and 
integrating backward in ^. 

The stable and unstable manifolds are best visualized by keeping track of their 
intersection with a Poincare (surface of) section V, a 3-dimensional surface in our 
4-dimensional phase space. Here we will be interested in symmetric homoclinic 
connections, which always intersect the symmetric section iS. It will therefore be 
important to choose a Poincare section that contains S, for instance ^ = 0. 

Figure [8] shows the first four intersections of our numerical approximation to the 
stable and unstable manifolds of Qq with V, for three different values of ui with fixed 
V = 0.8. For frequency value to = 0.11 (panel [8^) the stable and unstable manifolds 
intersect transversely and two homoclinic orbits (grey solitons) exist, marked by the 
points of intersection Ii (belonging to the homoclinic orbit in figure [2jd)- [2jf)) and 
I2 (belonging to the orbit in figure |8jd)). At w ~ 0.10545, Ii and I2 coalesce into a 
single solution and the manifolds become tangent. For uj < 0.10545 the unstable and 
stable manifolds do not intersect and the family of homoclinic orbits shown in figure 
|2] ceases to exist for V = 0.8. Note however that different families of solitary waves 
do exist for cj < 0.10545, V — 0.8; in figure [S] we only plot a part of the stable and 
unstable manifolds and therefore more intersections can still take place. This becomes 
apparent when comparing, for instance, figures [2] and [4] showing that a family of grey 
solitons can extent beyond the range of existence of another family. 




Figure 6. Some examples of symmetric multi-iiump grey solitons. 




Figure 7. Some examples of different types of solutions: (a) homoclinic 
asymmetric, (b) heteroclinic asymmetric and (c) heteroclinic with 2-humps. Panel 
(d) shows the projection of orbit (a) and (b) on the (j> — a plane. 
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A similar analysis was carried out for the second type of grey solitons, figure 
|4] We computed the stable and unstable manifolds for parameter values close to 
the the boundary exhibited in figure |4][a) and we did not observed a tangency 
close to the homoclinic orbit, confirming that the existence boundary is a numerical 
artifact. Therefore, the Poincare analysis is a useful tool to distinguish real boundaries 
from numerical artifacts. Note also that visualization of the stable and ustable 
manifolds of the fixed points on a Poincare surface of section is a very effective way of 
generating initial guesses for the numerical computation of homoclinic and hetcroclinic 
connections through the spectral method of section |4] 

o> = 0.11 V = 0.8 

(0 = 0.10545 V = 0.8 

(a) '-^ (b) 




Figure 8. Panels (a)-(c) show the intersection of the stable (black) and the 
unstable (grey) manifolds with the Poincare section (p = 0. Panel (d) is the grey 
soliton corresponding to the point labelled l2- 



6. Stability of the solitary waves 

Even though an exhaustive analysis of the stability of the solitary waves is beyond 
the scope of the present work, we sketch out here some relevant aspects. As 
previously mentioned, the electromagnetic circularly polarized wave is susceptible 
to the relativistic Raman and the modulational instabilities that would ultimately 
destroy the solitary wave. However, since the growth rates of these instabilities are 
controlled by the amplitude of the electromagnetic wave (oq in our dynamical system) 
and the plasma density (related to the parameter ui), certain parameter values could 
allow long distance propagation of the solitons. 

To test the stability of the solitary waves we show a couple of 1-dimensional 



simulations with the PIC code Calder 43 . This method does not prove stability 
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of the solutions in a strict mathematical sense but it provides an insight into the 
dynamics. We took a computational domain equal to 60 x 27r-\/l — V^/uiV (in c/ujpe 
units) that is large enough to assume periodic boundary conditions and we used one 
million of cells with ten particles per cell. The code was initialized with a grey solitary 
wave of the type presented in figure |4j 

Figure |9] shows the evolution of the electric field component of a grey solitary 
wave with V = 0.95 and uj — 0.95. For this high value of the frequency the ions are 
almost immobile and the asymptotic value of the amplitude of the wave is oq ^ 0.48. 
It propagates undistorted during a few tens of i^j^e^, until the part of the solution 
with vanishing E,j; (corresponding to the infinitely long circularly polarized wave) 
becomes unstable due to the Raman instability. However, a solitary wave with a 
lower frequency value lo = 0.5 (oq = 3.88) presents different dynamics (see figure 10). 
For this second case the solitary wave develops an instability at its trailing edge and it 
radiates part of its energy away. A similar instability appears in multi-hump solutions 
with VBC 44 . These examples reveal that, without collision and for a cold plasma 
model, our solutions present an unstable character. However, a dissipation mechanism 
could reduce or even supppress the Raman instability, thus allowing the propagation 
of the solitary waves for longer distances. 



60 

Time = 23.85 01 



(IJ) 









X (c/ta J 
Time = 27.9 k 



Figure 9. PIC simulation initialized with a grey soliton of the type shown in 
figure[4]with V = lu = 0.95. 



7. Conclusions 



Solitary waves excited by the interaction of a high-intensity laser with a plasma have 
been observed in laboratory experiments [6 -11 and particle-in-cell (PIC) simulations 



mm 



Although many theoretical works have been carried out on these structures 
(see 21 for a review), the parametric domain where the fixed point that controls the 
NVBC is a saddle-focus was unexplored. In our study of this regime we were able to 
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Figure 10. PIC simulation initialized with a grey soliton of the type shown in 
figure[4]with V = 0.95 and oj = 0.5. 



exhibit new ranges of solutions including grey and dark solitons. 

We recall that, for VBC, the solutions are organized in the uj — V plane on a set of 



infinitely many branches 17 . Each branch is characterized by the number of humps 
(or number of nodes) of the vector potential and it ends at a certain point due to the 
wavebreaking of the soliton. The potential is always a symmetric function whereas 
the vector potential can be either symmetric or antisymmetric. 

For NVBC we have shown that there is a continuum of solutions in the uj — V 
plane. Grey solitary waves with a symmetric potential and vector potential and dark 
solitary waves with symmetric potential and antisymmetric vector potential can be 
found for a wide range of parameters. This is a natural extension of the VBC case 
with the number of nodes being even or odd, respectively. Multi-hump solutions are 
also possible for any value of parameters for which a single-hump orbit exists. Further, 
asymmetric single and multi-hump solutions exist, consistent with symmetry breaking 
in a conservative, time-reversible system. 

In addition to being an important channel of laser-pulse energy transformation, 
solitary waves have been propossed as interesting candidates for photon and particle 



acceleration schemes 13 17 . PIC simulations with an overcritical amplitude soliton 
showed electron acceleration during the nonlinear wavebreaking [15] whereas ion 
acceleration has been detected during the postsoliton expansion pF In 17 , the 
authors reported wavebreaking of the solitary wave at the critical velocity determined 
by the end of the branch of the solutions in the uj — V plane. It was estimated that 
ions could reach an energy value of the order of 70 MeV. Similarly, in the present work 
we have presented domains of existence in the uj — V plane and we have shown that 
coalescence of solitary waves (in parameter space) leads to disappeance of families 
of solutions. Coalescence was visualized by keeping track of the stable and unstable 
manifolds of a fixed point by means of a Poincare surface of section. 
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We point out that the stabiHty of sohtons with VBC has been studied in the 
past [4| [T8||20{|23}|28] . In particular, 1-dimensional numerical fluid simulations with 
immobile ions showed that single hump solutions are stable whereas the multi-hump 
solutions suffer the Raman instability 27 44 . On the other hand, 2-dimensional 



simulations revealed that all solutions are unstable and the tranverse dynamics always 
dominates the longitudinal one 



28 



Our PIC simulations initialized with a grey 
solitary wave showed that some of them could propagate undistorted during a few 
tens of w~e^, just before the circularly polarized wave suffers the Raman instability. 
However, other grey waves radiate a portion of their energy from the trailing edge, 
similarly to the multi-hump solutions with VBC [44]. Since these are just a few 
examples, a complete stability analysis would requiry the study of other types of 
solutions (dark waves, asymmetric, multi- humps etc) in the whole cj — V plane. Adding 
warm plasma effects or a collision term would be relevant too. 

Besides existence and stability, the question about how to excite solitary waves 
with NVBC remains open. However, soliton-like electromagnetic modes with VBC 

Our preliminary PIC 



46 



have been observed during laser plasma interaction 
simulations on solitary wave excitation show that the interaction between a solitary 
wave with VBC and a long circularly polarized laser pulse can produce a solitary wave 
with NVBC. These simulations will be presented in a future work. 
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Appendix. The existence and stability of the fixed points 

Let us write (5a)-(56) as dx/dS, = fix) with x = [(fi (j) a d]. Fixed points are given by 



= a = and the solutions of 



V[- + - \ -co- 
re r,. 



V 



a = 
= 



whereas the stability depends on the eigenvalue of the Jacobian matrix 



[ 



V 







1 




-uj^ + v(^ + ^)+V{l-V^)a 



r'- 



The four eigenvalues of J at the fixed points Qq can be written as A^. 



with 

A(^,ao) 



> 



2T/2 



Ti 



{VTeV.Y 



(A.l) 
(A.2) 

° \ 



1 

V 



{A.3) 



(A.4) 
(A.5) 
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Fixed point Qq (or Qq ) is a saddle-focus if A < 0, a saddle-center if a/A > S > 
and a center if < \/A < S. The conditions a/A = S and A = yield the velocities 



(10 1 and (111 respectively. 



The eigenvalue of the Jacobian matrix at Qi for V — Vg are 



Al 2 = ±I^A 



/ (rer,-i)(re~r02 
2(re-r,)2 + (r, + ere)2a2 



(A.6) 



A., 



±i(l + e)= 



[(r, + 6,re)2-(i-v;2)(i + ,y 



13/2 



(A.7) 



and clearly the fixed point Qi is a saddle-center (note that T^jTi > 1). 

On the other hand, the existence and stability analysis of Q2 requires some 



auxiliar operations. From (A.l) and assuming — Fg <<pf < F^/e, one gets 

2 i^l-i^l 



and by substituting in (A.2) 



G{4>f) 



1 



F F 



V 



fpetpi V Ipt 



.2^2_(l_y2)(l_£2) 



= 



(A., 



(A.9) 



where we introduced the subscript / to denote that we are dealing with a fixed point. 
Note that restriction ([7| together with 
must lie on the intervals 



> in (A.8) show that solutions of (A.9 1 



F — F- 



<(l)f< 



2e(F,-|-eFe) ~ 



(A.IO) 



The solutions of (A.9) can be discussed taking into account some properties of 
the function G{(j)f) and its derivative: 



dG 



V 



2^1 



2^1 - (1 - V2)(l - £2) 



i'elpi 



2^^ 



-(1 



- V2)(l 



5!L 

.2V 



{■ll>i - ti)e) 



In particular we are interested in the zeros of this derivative. The factor inside the 
squared root vanishes &t (jjf — (j>max2/V'^ > 4>max2, that it is outside the physical 
domain. On the other hand, the roots of the term inside the braces are given by the 
zeros of the cubic equation 



8(r, 



'§^e^<t>]-^{Z + V'^)e^^'} + i{V,-eT^)(l + V^)e4 



-{Ti-eTefv'^ + e{l-V'^)reTi = (A.12) 



For discussing the solutions of (A.9) within the domain 
we first note that G(0) 



< <t>l < 



^max2 ! 



0, corresponding with the fixed points Qq. One also has 
the asymptotic behaviours G — > —00 as 0/ — (f>max2 and G{(j)min2) > (< 0) for 
velocities less (greater) than v (see (13)). The analysis of its derivative shows that 



G has or 1 extreme for V < v and therefore (A.9) has one solution in this regime 
(Qq ). On the other hand, for > w it always has one extreme, at say 0^ with 
G(0p > 0. Therefore, if G{(/)*j:) > 0, (A.9) has two solutions and, in addition to QJ, 



there is another fixed point that we call Q2- One also checks that Q2 and Qg lie in 

Hence, heteroclinic connections among Q2 and Qq 



different manifolds given by (|6| 
are not possible. 



(A.ll) 
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